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NOTATION 


The symbols used in this document and their definitions are listed below for 
convenience. 


Roman Symbols 

a. . . speed of sound 

Cp . . . specific heat at constant pressure 
Cy . . . specific heat at constant volume 
e . . . internal energy 

1 . . . z index of numerical solution 
j ... r index of numerical solution 

k . . . 0 index of numerical solution or thermal conductivity 

1 . . . turbulence model damping function 
n . . . outward unit normal vector 

p . . . pressure 

r . . . radius or radial coordinate 
t . . . time 
v . . . velocity 
z . . . axial coordinate 
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A. . . surface area 

A+ . . . turbulence model parameter 
Cep • • • turbulence model parameter 
Cki e b • • • turbulence model parameter 
C wa f te • • • turbulence model parameter 

CFL... Courant-Freidrichs-Levy number (At/ '^t maXr $ t a ble) 

D ... dissipation flux vector, turbulent damping parameter, or diameter 
F . . . flux vector in z direction or turbulence model function 
G . . . flux vector in r direction 
H . . . flux vector in 0 direction 
H i . . . total enthalpy 

K . . . source term flux vector or turbulence model parameter 

L . . . length 

M . . . Mach number 

Pr . . . Prandtl number 

Pr turbulent ■ ■ • turbulent Prandtl number 

Q . . . vector of dependent variables 

R. . . gas constant or residual 

S . . . blade row correlations or pertaining to surface area normal 
T ... temperature or torque 
V . . . volume 

Greek Symbols 


viii 


a . . . time-stepping factor 



. . . modified second-order damping coefficient 
. . . modified fourth-order damping coefficient 
p . . . density 

k? . . . second-order damping coefficient 
. . . fourth-order damping coefficient 

7 . . . specific heat ratio 

6 . . . spatial second-order central difference operator 
\ . . . blockage factor 

c\ 

X v . . . second coefficient of viscosity (= — j/i) 
p . . . coefficient of viscosity 

T] . . . dimensionless wall normal coordinate (= y\J ) 
v . . . damping factor 
u> . . . anglular velocity or vorticity 
A . . . increment of change 

Special Symbols 

V . . . spatial vector gradient operator 

A . . . spatial forward difference operator 

V . . . spatial backward difference operator 

Superscripts 


[ ] . . . averaged variable 

[ ] . . . dimensional variable 

. . . implicitly smoothed variable 
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[~ *] . . . vector variable 
[ ]*... intermediate variable 

[ ] n . . . time step index of variable 

Subscripts 

[ \ef fective • • • effective flow value 
[ ]» j fe • ’ • point index of variable 
[ )i nv . . . inviscid component 
[ haminar • • • laminar flow value 
[ }max • • • maximum value 
( 1mm v • minimum value 
[ ]p . . . related to pressure 
[ )t • • • total quantity 

[] z . . . derivative or value with respect to z 
[ ] r . . . derivative or value with respect to r 
[ . . , derivative or value with respect to 8 

[ turbulent • • ■ turbulent flow value 
[ ]qo . . . freestream value 
[ ] re y . . . reference value 
[ • • • Klebanoff intermittency factor 

[ \vis • • • v i sc °us component 
[ ] wa } te . . ■ turbulent flow wake parameter 
[ ]2 . . . second-order value 
[ ]4 . . . fourth-order value 



1. SUMMARY 


i 


The purpose of this study is the development of a three-dimensional Euler/Navier- 
Stokes flow analysis for fan section/engine geometries containing multiple blade rows 
and multiple spanwise flow splitters. An existing procedure developed by Dr J. J. 
Adamczyk and associates and the NASA Lewis Research Center was modified to ac- 
cept multiple spanwise splitter geometries and simulate engine core conditions. The 
procedure was also modified to allow coarse parallelization of the solution algorithm. 
This document is a final report outlining the development and techniques used in the 
procedure. 

The numerical solution is based upon a finite volume technique with a four stage 
Runge-Kutta time marching procedure. Numerical dissipation is used to gain solu- 
tion stability but is reduced in viscous dominated flow regions. Local time stepping 
and implicit residual smoothing are used to increase the rate of convergence. Mul- 
tiple blade row solutions are based upon the average-passage system of equations. 
The numerical solutions are performed on an H-type grid system, with meshes being 
generated by the system (TIGG3D) developed earlier under this contract. The grid 
generation scheme meets the average-passage requirement of maintaining a common 
axisymmetric mesh for each blade row grid. 


1 



I The analysis was run on several geometry configurations ranging from one to five 
blade rows and from one to four radial flow splitters. Pure internal flow solutions were 
obtained as well as solutions with flow about the cowl/nacelle and various engine core 
flow conditions. The efficiency of the solution procedure was shown to be the same 
as the original analysis. 
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2. INTRODUCTION 


This document contains the Final Report for the ADPAC-APES (Advanced 
Ducted Propfan Analysis Codes- Average Passage Engine Simulation) program devel- 
oped by the Allison Gas Turbine Division of the General Motors Corporation under 
Task IV of NASA Contract NAS3-25270. The objective of this task is development of 
a three-dimensional flow analysis tool for advanced fan section and turbofan engine 
geometries such as the NASA/GE Energy Efficient Engine seen in Fig. 2.1 . The tool 
is able to compute steady flow' solutions about geometries with any number of blade 
rows and axisymmetric radial flow 7 splitters. The tool computes the flow through 
the fan and optionally about the fan cowl and engine nacelle, both upstream and 
downstream of the engine. When the domain is extended in this manner, engine 
performance can be determined entirely by the analysis tool. Effects of engine core 
flow can also be simulated. 

Details of the flow solution algorithm are covered in Chapter three of this docu- 
ment. Chapter four presents solution results for various geometries and comparisons 
to experimental data. A summary of the conclusions for this study is given in chapter 
five. 

This flow' analysis tool was developed from a code entitled VSTAGE which was 
developed by John J. Adamczyk of the NASA Lewis Research Center [1] . The user is 
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referred to the documentation for that code for additional information on the solution 
procedure. 
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Figure 2.1: NASA/GE Energy Efficient Engine cross section 
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3. EULER/NAVIER-STOKES NUMERICAL ALGORITHM 


This chapter describes the algorithm for the flow solver and outlines a solu- 
tion procedure for multiple blade row calculations. As stated earlier, the solver was 
primarily developed by J. J. Adamczyk [1] [2] which was based upon a procedure 
originated by Jameson [3]. The definitions of the pertinent variables used in this 
chapter may be found in the Nomenclature. 


3.1 Nondimensionalization 


The variables in the numerical solution are nondimensionalized by reference val- 
ues as follows: 


z — 


J ref 


r = 


f v z v r vq 

— > v * - z . — ■> v r = z — :> v e = — 


P = 


Pref 


P = 


■'ref 


P 


T = 


Pref 

T 


C P = 


v ref 
c p 


R 


P = 


ref 

P 


c v = 


u> = 


v ref 

cy 

^ref 
< * >J Wef 


k = 


v ref 

k 


*ref 


1 ref Pre f v ref 

The reference quantities are defined as follows: 

L re j The maximum diameter of the blade represented in the grid 
p re f The reference (or freestream) relative total pressure 
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p re j The reference (or freestream) relative total density 

t; f The reference (or freestream) velocity determined from the 
rej 

relative total conditions v re j = \JlP re f I P T ef 
p re f The reference (or freestream) viscosity 
fc re y The reference (or freestream) thermal conductivity 
R re j The reference (or freestream) gas constant 
T re j The reference (or freestream) temperature 


3.2 Governing Equations/Discretization 

The numerical solution procedure is based on the strong conservation law form of 
the Navier-Stokes equations expressed in a cylindrical coordinate system. The Euler 
equations may be derived as a subset of the Navier-Stokes equations by neglecting vis- 
cous dissipation and thermal conductivity terms (i.e. - /i and k = 0). A derivation of 
the average- pas sage equation system can be found in [4]. For a multi-blade row envi- 
ronment, this equation system is obtained by filtering the Navier-Stokes equations in 
time and space to remove all information except the time average flowfield pertaining 
to a specific blade row'. For the particular blade row, integration of these equations 
equations over a rotating finite control volume produces the following equations: 


/ j t (\Q)dV + L inv (\Q) = / A SdV + / A KdV + L vis (XQ) (3.2) 
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The vector of dependent variables Q is defined as: 

P 

pvz 

pv r (3.3) 

P v 6 

PH- 

where the velocity components Vz,vr , and vq are the absolute velocity components 
in the axial, radial, and circumferential directions of the coordinate system for the 
fan section (see Fig. 3.1), respectively. The term A represents the neighboring blade 
row blockage factor This factor has a value between zero and one, with a value of one 
indicating zero thickness for neighboring blades. 

The term A S contains the body forces, energy sources, and momenta correlations 
associated with the neighboring blade rows. The terms L^ nv and L v j s represent the 
cell face mass, momentum, and energy flux evaluations for the inviscid, and viscous 
components, respectively. These terms are defined as: 

= j iA [a F inv iA z + A O inv iA r + A (H inx> - r wQ)dA e ] (3.4) 

and: 

£« S (A<?) = J iA [A Fvi s iA z + A G vU dAr + A (H vis - ru,Q)dA„] (3.5) 

The individual flux functions are defined as: 
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Cartesian Reference 
Coordinate System 


Figure 3.1: Fan section/engine analysis computational domain 
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F = F(Q), 0 =G(<5), H = H(Q) 


Fv = F v (Q), 


Gv = G„(l?), Hv = H V (Q) 


( 3 . 6 ) 


( 3 . 7 ) 


( 3 . 8 ) 


The flux variables F,G, and H are determined at each grid cell interface by deter- 
mining the average (Q) of the cell-centered dependent variables from the individual 
finite volumes adjoining the interface. 

Finally, the cylindrical coordinate system source term is: 


K = 


0 

0 


T ee 


o 

o 


( 3 . 9 ) 


It should be noted that in the numerical algorithm, the radius used in the cylindrical 
source term K is carefully formulated to guarantee numerical conservation for the 
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radial momentum equation. That is, for a uniform stagnant flow, the radius in the 
radial momentum equation is chosen such that both sides of the radial momentum 
equation are equal. This ensures that small geometric errors do not corrupt the 
conservative nature of the numerical scheme. 

The total energy function, e^, is defined as: 

= + ^ + ^ + ^ (310) 
The total enthalpy, H, is related to the total energy by: 

H = e t + V - (3.11) 

P 

The viscous stress terms may be expressed as: 


Tzz — 2/i 


T Z r = P 


T zd = 2 P 


' dv£\ 
K dz) 

' dv r \ 

1 dvr 
r 89 

dv r 


+ A v V-V, 


+ 


'dvz 
, dr 


+ 


(t*)] 


T rr = 2p I 1 + A-y V • V, 


— |(S) *(£)-(?)] ■ 

T ee = 2 p 


1 8vq 8v r 
r 89 ^ r 


+ A v V • V, 


8T 


qz = v z Tzz + v r T Z r + V Q T Z 9 + 

8T 

q r = VzTtz + VrTrr + v Q T r9 + 


(3.12) 

(3.13) 

(3.14) 

(3.15) 

(3.16) 

(3.17) 

(3.18) 

(3.19) 
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(3.20) 


. dT 

QQ = v ztq z + v r r er + VQT 6e + k — , 
where /x is the first coefficient of viscosity, is the second coefficient of viscosity, 


and: 


V V — — 1<9u fl 

dz + dr + r 39 r 


(3.21) 


The remaining viscous stress terms are defined through the identities: 


Trz = T Z r, 


( 3 . 22 ) 


T 6r T rd ’ 


(3.23) 


T 6z ~ T z6 ’ 


(3.24) 


This integral form of the governing equations is applied to a generalized finite volume 
in physical space as shown in Fig. 3.2. The cell surface areas dAz , dAr , and dAg are 
calculated using the cross product of the diagonals of a cell face, and the cell volume 
is determined by a procedure outlined by Hung and Kordulla [5] for generalized 
nonorthogonal cells. 


3.3 Fluid Properties 


The fluid is assumed to be air acting as a perfect gas, thus the ideal gas equation 
of state has been used. Fluid properties such as specific heats, specific heat ratio, and 
Prandtl number are assumed to be constant. The fluid viscosity is derived from the 
Sutherland formula: 


/* = Cl 


3 

t + c 2 


(3.25) 
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• Grid Point 

Grid Line 

Hidden Grid Line 



Figure 3.2: Three-dimensional finite volume cell 
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The so-called second coefficient of viscosity \ v is fixed according to: 



(3.26) 


The thermal conductivity is determined from the viscosity and the definition of the 
Prandtl number as: 



(3.27) 


3.4 Artificial Dissipation 

The discretized system of equations has unstable properties and can exhibit odd- 
even point decoupling, especially near regions where high gradients of the flow quan- 
tities Q exist (e.g., shocks). To suppress these instabilities, an artificial dissipation 
operator ( D ) is added to the numerical scheme. Jameson [3] demonstrated that a 
dissipative system combining second- and fourth-difference smoothing terms can ef- 
fectively eliminate undesirable numerical oscillations without destroying the accuracy 
of the solution. 

The dissipation operator for the first index is shown below: 


DZ ^ d i+y,j,k d i-y,j,k 


(3.28) 


d. i 


where: 


V j 


(3.29) 




(3.30) 


\ 


15 



_ IPi + lJ,fc- 2 P»J,fc+P»-lJ > fcl 

Vl ' >J ' k \Pi+l,j,k + 2 Pi,j,k + Pt-1 J,fcl 

Typical values for the second and fourth difference damping constants are: 


(3.31) 

(3.32) 




k 4 = 


64 


(3.33) 


The term represents a one-dimensional equivalent of the maximum allowable 
time step in the given coordinate direction. The use of this factor introduces an 
eigenvalue scaling into the dissipation operator which minimizes the added dissipation 
in coordinate directions which do not limit the stability of the algorithm. 

The scheme presented above is stable for all time steps satisfying the C F L - related 
time step limitation 

CFL < 2V2 (3.34) 


The damping scheme described above may be applied directly for inviscid flow 
calculations, but must be modified slightly for viscous flow calculations. In regions of 
the flowfield where viscous dissipation increases, the artificial dissipation should be 
reduced. In order to produce this effect, a Mach number scaling is employed and the 
modified second and fourth difference coefficients are shown below. 


M. 1 

o . . , "^iJik 

{€2) i+y,k = K max{u i+hhk ’ Jj ’!) 


(3.35) 


M. 1 . , 


(3.36) 


where M is the mid passage or free stream Mach number. 
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The complete dissipation operator D j j ^ is constructed as the sum of the dissi- 
pation operators in each of the respective coordinate directions as: 




(3.37) 


3.5 Time Integration 


The time-stepping scheme used to advance the discretized equations is a four- 
stage Runge-Kutta integration. The solution proceeds as: 

Q 1 = Q n -a 1 At[L(Q n )+D(Q n )l 

Q 2 = Q n -a 2 At[L(Q 1 ) + D(Q n )}, 

Q 2 = Q n -a 3 At[L(Q 2 ) + D(Q n )), 

Q A = Q n -a i At[L(Q 2 ) + D(Q n )], 


<? n+1 = 

where: 

1 1 1 1 

a l = g> a 2 = 4> a 3 = 2’ a 4 = 1 


and: 


£«?) = U«)-I.i,(«) 


(3.38) 


(3.39) 


(3.40) 


Linear stability analysis indicates that this scheme is stable for all time incre- 
ments At which satisfy the stability criteria CFL < 2\/2. The CFL number may be 
defined in a one-dimensional manner as: 


CFL = 


At 
u -fa 

~kr 


(3-41) 
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This factor is calculated for each coordinate direction, and then geometrically av- 
eraged to obtain the maximum allowable time increment for a given computational 
cell. 

The acceleration technique known as local time stepping is used to enhance 
convergence to the steady-state solution. Local time stepping utilizes the maximum 
allowable time increment at each finite volume cell during the course of the solution. 
If a truly unsteady flow calculation is desired, a uniform value of the time step A t 
must be used at every cell to maintain the time-accuracy of the solution. 

3.6 Implicit Residual Smoothing 

Implicit residual smoothing is used to extend the stability limit of the algorithm 
and increase the rate of convergence to a steady state solution. Residual smooth- 
ing attempts to accelerate the propagation of changes in the dependent variables by 
filtering the residuals of the calculation (which may also be interpreted as the local 
time derivative of the computational solution) at each time step. By enhancing the 
transfer of information between grid points, calculation time steps much larger than 
the stability-limited values may be utilized. Residual smoothing was originally intro- 
duced by Lerat (see e.g. Hollanders, et al. [6]) for use with the Lax-Wendroff scheme 
and later applied to Runge-Kutta schemes by Jameson and Baker [7] as a technique 
to accelerate convergence for steady-state calculations. 

Since the time rate of change of the dependent variables dQ/dt is in essence 
controlled by a residual operator R(Q) = L(Q) — D(Q) ) it would follow that any 
measure which accelerates the propagation of changes in the residual throughout 
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the domain would ultimately enhance convergence. The implicit residual smoothing 
operator used in this study can be written as: 

(1 - e z S Z z)( 1 - «r£rr)(l - = ^i,j,k (3.42) 

where the differencing operator 8 is expressed as: 


8 zz R^ j k - R-i + lj,k 2 R i,j,k + ^i—l,j,k 


(3.43) 


A value of e zz = trr = £QQ = 2 is typically used. 

The reduction is applied along each coordinate direction separately as: 


H,j,k = (1 ~ lR ij,k 

(3.44) 

!**,, = (! -..M-Xvk 

(3.45) 

= a - 

(3.46) 

Kit = KXk 

(3.47) 


where each of the first three steps above requires the inversion of a scalar tridiagonal 
matrix. The residual smoothing operator is applied to the first and third stage of the 
four-stage Runge-Kutta algorithm. The time-marching scheme then becomes: 


Q l =Q n -a 1 R(Q n ) 


Q 2 = Q n -a 2 R(Q i) 


Qz = Q n - <* 3 R{Q 2 ) 


Q 4 = Q n -cx A R(Q z ) 

Q n+1 = Qa 


(3.48) 
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The implicit residual smoothing operator applied in this context allows a time step 
greater than the unsmoothed stability limited step. For example, the new time step 
limit for the axial coordinate must satisfy: 


1 


( ^_^^_svnooth \ 
V CFL ) 


(3.49) 


Thus with the unsmoothed stability criteria ( CFL < 2\/2), a sample of new limits is 
listed below. 


e z = 1.0 

CF^mootb £ «- 32 

O 

cs 

II 

b* 

CFL smooth < 8.48 

e z = 3.0 

CFL smoo th < 10.19 


3.7 Turbulence Model 

The effects of turbulence are accounted for with a relatively standard version of 
the Baldwin- Lomax [8] turbulence model. This model is computationally efficient, 
and has been successfully applied to a wide range of geometries and flow conditions. 

The effects of turbulence are introduced into the numerical scheme by utilizing 
the Boussinesq approximation , resulting in an effective calculation viscosity defined 
as: 

^ef f ective = P laminar f 1 turbulent (3.50) 

The simulation is therefore performed using an effective viscosity which combines the 
effects of the physical (laminar) viscosity and the effects of turbulence through the 
turbulence model and the turbulent viscosity /x turbulent • 
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(3.51) 


The Baldwin-Lomax model specifies that the turbulent viscosity be based on an 
inner and outer layer of the boundary layer flow region as: 

j ( Hurlmlent)inner > V ^ Vcrossover 

P turbulent = ) , x 

l \PturbulenV outer •> V > Vcrossover 

where y is the normal distance to the nearest wall, and ycrossover is the smallest 
value of y at which values from the inner and outer models are equal. The inner and 
outer model turbulent viscosities are defined as: 


(pturb^inner ~ P^ W 

(Pturb^outer = ^■^' c VP F wake F kleby 
Here, the term l is the Van Driest damping factor 

l = ky(l-e(-y + / A+ )) 

u> is the vorticity magnitude, F wa ^ e is defined as: 

F wake = VmaxFmax 

where the quantities ymax, Fmax are determined from the function 

F{y) = 3/Ml 1 - e ^~ V 

The term is defined as 


(3.52) 

(3.53) 


(3.54) 


(3.55) 


y, 


p M 


(3.56) 


(3.57) 


\ M laminar 

The quantity Fmax is the maximum value of F(y) that occurs in a profile, and ymax 
is the value of y at which it occurs. The determination of Fmax an d ymax is perhaps 
the most difficult aspect of this model for three-dimensional flows. The profile of F(y) 
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versus y can have several local maximums, and it is often difficult to establish which 
values should be used. It has been found from numerical experimentation that the 
most reliable value of Fmax is taken as the maximum value of F(y ) between a y+ 
value of 350.0 and 1200.0. The function is the Klebanoff intermittency factor 

given by 

u>) = [i + 5 -«(^) 6 r 1 < 3 - 5 *) 

and the remainder of the terms are constants defined as: 


A + = 26, 

Cep = 1 . 6 , 

C kleb = °- 3 ’ 
k = 0.4, 

K = 0.0168 (3.59) 

In practice, the turbulent viscosity is limited such that it never exceeds 1000.0 times 
the laminar viscosity. 

The turbulent flow thermal conductivity term is also treated as the combination 
of a laminar and turbulent quantity as: 


^ef fective ~ ^ laminar ^ turbulent (3.60) 

For turbulent flows, the turbulent thermal conductivity k^ ur ^ u i en ^ is determined from 
a turbulent Prandtl number P r turbulent suc ^ ^at 


Pr 


turbulent 


^turbulent 


(3.61) 
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The turbulent Prandtl number is normally chosen to have a value of 0.9. 

In order to properly utilize this turbulence model, a fairly large number of grid 
cells must be present in the boundary layer flow region, and, perhaps of greater im- 
portance, the spacing of the first grid cell off of a wall should be small enough to 
accurately account for the inner “law of the wall” turbulent boundary layer profile 
region. This requires the first cell in the laminar sublayer with a y~^~ value typically 
around 5. Unfortunately, this constraint is typically not satisfied due to grid-induced 
problems or excessive computational costs, especially for multiple blade row calcu- 
lations. Practical applications of the Baldwin-Lomax model for three-dimensional 
viscous flow must be made with the limitations of the model in mind. The Baldwin- 
Lomax model was designed for the prediction of wall bounded turbulent shear layers, 
and is not likely to be well suited for flows with massive separations or large vortical 
structures. 


3.8 Boundary Conditions 


Inflow and exit boundary conditions are applied numerically using characteristic 
theory. A one-dimensional isentropic system of equations is utilized to derive the 
following characteristic equations at an axial inflow /outflow boundary: 


8C~ 

dt 


- (v z - a) 


dC~ 

dz 


= 0 , 


dC+ 

dt 


+ ( v z + a) 


<9C+ 

dz 


= 0 


(3.62) 

(3.63) 


where: 

C~ = VZ - ——t, C+ = V Z + — (3.64) 

7-1 7-1 


23 



These boundary condition equations are based upon an inflow inlet and an outflow 
exit. Aerodynamic conditions not satisfying these requirements (e.g., reverse flow) will 
cause spurious results or failure. In order to efficiently process boundary information 
in the numerical solution, phantom cells are located just outside the computational 
domain to permit the unmodified application of the interior point scheme at near 
boundary cells. 

For subsonic normal inflow, the total pressure, total temperature, and flow angle 
are specified. The upstream running invariant C~ is extrapolated to the inlet, and 
the equation of state and flow equations are used to determine the variables at the 
inlet boundary. 

At the exit, a static pressure is specified at the hub for internal flows, and at 
the outer boundary for external flows. The remaining pressures along the outflow 
boundary are calculated by integrating the radial momentum equation: 

?P = ^£ (3.65) 

dr r 

In this case, the downstream running invariant C+ is used to update the phantom 
cells at the exit boundary. Far-field boundaries also use this characteristic technique 
based on whether the local flow normal to the boundary passes into or out of the 
domain. 

For applications where core flow conditions are simulated (e.g., combustor , high 
pressure spool device) boundary conditions similar to those just discussed are em- 
ployed. The entrance to the core region is treated as a local exit of the domain, even 
though the region can be inside the computational domain. An example of this is 
seen in Fig. 3.3 . A specified hub static pressure, radial momentum equation, and 
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Figure 3.3: Computational domain highlighting core definition 

characteristic equation set the flow variables. Similarly, the exit to the core region is 
treated as an independent inflow of the domain. Specified core flow total pressure, to- 
tal temperature and flow angle are used with the characteristic equation to determine 
the flow quantities at the core exit boundary. 

All solid surfaces (hub, cowl, radial flow splitters, and airfoils) must satisfy flow 
tangency for inviscid flow: 


V -n = 0 


(3.66) 
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or no slip for viscous flows: 


v z = 0, v r = 0, vq = tu> (3.67) 

In both cases, we specify no flux through the boundary (an impermeable surface), 
and hence, only pressure is needed at the phantom cell. The pressure is determined 
by extrapolation. Solid surfaces are also assumed to be adiabatic, which implies that 
the normal temperature gradient is also zero. 

The calculation presumes that the flow solution is periodic with a period of 
one pitch (arclength between tangential extrema of the grid). Therefore, all cells at 
the tangential boundaries of the domain (and not defining a solid surface) take as 
their phantom cell flow variables the quantities from the cell volume at the opposite 
tangential bound. 


3.9 Solution Procedure 


A procedure for obtaining a numerical solution for multiple blade rows is de- 
scribed below. The single blade row case is in general a reduction of the multiple 
blade row case and is described later. Before executing the solution algorithm, nu- 
merical grids (one for each blade row) are required. These grids model the actual 
three-dimensional geometry of their particular blade row, and represent the rest of 
the domain as an axisymmetric duct. The average passage method requires all grids 
to have the same meridional representation (i.e., the same dimensional (z,r) coor- 
dinate lattice structure). More information on the required grids can be found in 
[9]. 
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Once the grids are obtained, the solution procedure is begun from a specified 
initial condition. This initial condition is a uniform flow, or is introduced from a 
previous solution. The multiple blade row solution is found using a nested iteration 
procedure with an inner and outer loop as seen in Fig. 3.4. Within the inner loop, 
the Runge-Kutta time integration is used and the Euler/Navier- Stokes equations are 
solved for a particular blade row (grid). As mentioned earlier, the average passage 
form of the equations are used and the neighboring blade row effects (blade forces, 
correlations) are modeled as steady parameters for the inner loop time integration. 
When each blade row in the domain has gone through this iteration and the blade 
row effects have been recalculated, one cycle or “flip” through the system (outer loop) 
is complete. Once during a flip, each blade row’s force terms are updated based upon 
the the axisymmetric average of that blade row’s flow field. When these updated 
effects are used depends upon user control of the solution procedure. For Fig. 3.4 the 
neighboring blade row effects are all from the previous flip. However, if the updated 
terms for the current flip are used, the solution procedure is represented in Fig. 3.5. 
For more information on solution procedure techniques see [10]. Solutions are deemed 

O 

converged when the average residual R has been reduced by a factor of 10 

For a single blade row case, the solution procedure is greatly simplified. With 
no neighboring blade rows and their effects to calculate, there is no outer loop. The 
Runge-Kutta time integration for the equation system of the single blade row is 
executed until the convergence criterion is met. 
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Figure 3.4: Multiple blade row solution procedure 
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BLADE ROW # 2 ROWS 


Figure 3.5: Solution procedure using updated blade row effects 
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4. RESULTS 


In the following sections, solutions from the computational procedure described 
in the previous chapter are presented. An Euler solution for the fan section of the 
GMA 3007 will be presented and compared to experimental data. Also, an Euler 
solution for the GE/NASA Energy Efficient Engine (E cubed) fan section will be 
shown. Following this, a flat plate test case showing the boundary layer characteristics 
predicted by the Navier-Stokes algorithm will be presented. Finally, comparison of 
a viscous solution for the GMA 3007 fan section to experimental data will conclude 
the chapter. 


4.1 GMA 3007 Fan Section - Euler Analysis 


The initial verification of the original Euler flow solver was presented in [1] and 
the solver has essentially remained the same since then. To test the multiple flow 
splitter capability, calculations were done on the geometry of the GMA 3007 fan 
section. The simulation of the fan section included the rotor, the core duct guide 
vane, and the bypass duct guide vane. One radial flow splitter (the core/bypass duct 
splitter) is part of the geometry, and this splits the exit of the computational domain. 
All these features can be seen in Fig. 4.1 which shows an axisymmetric plane of the 
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Table 4.1: Grid parameters for the GMA 3007 geometry 


Blade Row 
Grid 

Streamwise 

Points 

Spanwise 

Points 

Tangential 

Points 

Blade Pts 
LE to TE 

Blade Pts 
Hub to Tip 

Fan Rotor 

101 

29 

15 

21 

29 

Core Duct 
Vane 

101 

29 

11 

13 

17 

Bypass Duct 
Vane 

101 

29 

11 

17 

17 


computational grid. The full grid consists of 101 streamwise points, 29 spanwise 
points, and 15 tangential points. Some sections of the full grid for the rotor 
are shown in Fig. 4.2. Both the core vane and bypass vane grids are (101X29X11). 
Table 4.1 shows pertinent statistics for the grids used. The grid does not model the 
rotor tip clearance region. 

Test data for this geometry was available and Fig. 4.3 displays the locations 
where data was taken on a schematic of the fan section. The data consisted of total 
pressure and temperature measurements taken at the leading edges of both vane rows 
and at radial rakes downstream of the vanes. 

The fan section was simulated with the rotor at design speed. Converged solu- 
tions were obtained for different total pressure ratios by setting the exit static pres- 
sure condition. The performance characteristic (mass-averaged total pressure ratio 
vs. mass flow) of these solutions is compared to the experimental values in Fig. 4.4. 
The Euler solutions consistently over flow at this rotor speed and this is primarily 
caused by the lack of boundary layer blockage on the blade row surfaces. At the 
highest pressure ratio condition (near stall) the difference in mass flow is greatest and 
again this is not atypical of inviscid analyses. The shock-boundary layer interaction 
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Figure 4.1: Axisymmetric plane of GMA 3007 Euler grid 





Figure 4.2: 3-Dimensional grid for GMA 3007 rotor grid 
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which usually has a significant effect on the performance of a highly loaded blade row 
is not modeled in these solutions. 

Mass-averaged radial profiles of total pressure ratio are shown for two operating 
conditions of the fan section in Fig. 4.5. The pressure ratios are compared at the 
Rotor Exit Plane which is seen in Fig. 4.3. No data was actually taken at this plane 
and the test data shown in Fig. 4.5 has been back calculated from the vane leading 
edge locations based on predicted streamlines. Comparing results at the true data 
locations shows the same characteristics seen in Fig. 4.5. For both conditions, the 
Euler solution produces more pressure ratio in the lower portion of the span and less 
pressure ratio in the upper portion, with this discrepancy decreasing at the near stall 
condition. However, it should be noted that several Navier-Stokes analyses have been 
run on this fan section and all show more pressure rise near the hub than that seen in 
the test data. A potential explanation is that there is an error in the test data. The 
discrepancy in the upper portion of the span can be related to low grid resolution, the 
lack of a tip clearance model, and the inability of the inviscid solution to accurately 
model the boundary layer effected shock behavior in the blade passage. 

Figure 4.6 compares predicted and measured efficiency profiles at an operating 
condition near the design point. While the Euler solution’s mass-averaged efficiency 
is similar to the test data, the radial profile is significantly different. Part of this 
difference is likely due to the inviscid analysis’s lack of a model for end wall and tip 
clearance losses. 

Static pressure contours for the rotor flowfield are shown in Fig. 4.7. Pressure is 
displayed on the solid surfaces of the rotor grid and the blade surface pressure for the 
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Figure 4.4: GMA 3007 design speed characteristic 
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Figure 4.5: GMA 3007 Euler solution: total pressure ratio profile 









Figure 4.7: GMA 3007 Euler solution: color contours of static pressure 
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two vane rows are included as well. A portion of the core/bypass duct splitter has 
been removed to display more of the core vane row. The passage shock’s termination 
on the rotor pressure surface is seen on the upper blades. The passage shock was in 
the aft portion of the blade row which is typical of inviscid simulations. The influence 
of the splitter on the rotor performance was not significant and that will be discussed 
further in a later section. 

4.2 Energy Efficient Engine Fan Section 

The flowfield about the GE/NASA Energy Efficient Engine (E cubed) fan section 
was also simulated to test the Euler version of APES. The geometry included the fan 
rotor, a booster stage, a core duct guide vane, and a bypass duct guide vane. The 
geometry can be seen in Fig. 4.8, which is an axisymmetric plane of the computational 
grid. There are three spanwise flow splitters. A core/bypass splitter breaks up the 
exit of the domain into two regions. There is an island splitter over the booster stage, 
and there is also a part span shroud which splits the flow on the fan rotor. The grid 
for each blade row has 185 streamwise points and 45 spanwise points. Table 4.2 shows 
pertinent statistics for the grids used. 

This fan section geometry was tested and the reader is referred to [11] and [12] for 
details of the geometry and experimental configuration. The primary aerodynamic 
data taken was a performance map for the two air streams divided by the core/bypass 
duct splitter. For an operating condition near the design point, spanwise and tangen- 
tial measurements were available as well as static pressures on the endwall surfaces. 
Figure 4.9 shows a schematic of the experimental measurement locations. 
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Figure 4.8: Axisymmetric plane of E cubed Euler grid 



Table 4.2: Grid parameters for the E cubed geometry 


Blade Row 
Grid 

Streamwise 

Points 

Span wise 
Points 

Tangential 

Points 

Blade Pts 
LE to TE 

Blade Pts 
Hub to Tip 

Fan Rotor 

185 

45 

15 

29 

45 

Booster Stage 
Vane 

185 

45 

11 

15 

19 

Booster Stage 
Rotor 

185 

45 

11 

15 

19 

Core Duct 
Vane 

185 

45 

11 

13 

11 

Bypass Duct 
Vane 

185 

45 

15 

33 

33 



Figure 4.9: E cubed experimental measurement location schematic 


A Euler solution was obtained at a condition near the design point and the 
mass-averaged performance for the flowfield is shown on the experimental data maps 
of Figures 4.10 and 4.11. The bypass duct value showed reasonable agreement with 
the experiment, with the solution overflowing by 2%. The core duct data showed 
significantly less mass flow than the experimental value. This discrepancy is due in 
part to the core vane row in the solution not representing the actual geometry. The 
core duct vane geometry was not available so an approximate model for the vane row 
was used. 

Spanwise measurements were taken downstream of the bypass duct vane and 
this data is compared to the solution in Figures 4.12 through 4.14. Total pressure 
ratio is compared in Fig. 4.12 and the total temperature in 4.13. The range bars 
in the figures represent the variation in the stagnation quantity as recorded by the 
tangentially spaced arc rakes. The low values in the range are the wake measurements 
and the circle symbols represent mass averages of the arc rake data. Through much of 
the span, the comparison is reasonably good. However, at both endwalls deviations 
are apparent. The tip clearance flow not being modeled might be the cause of the 
discrepancy near the case. The missed performance in the core stream accounts for 
the differences seen near the splitter surface. The spanwise profile of efficiency is 
shown in Fig. 4.14. Again, the solution deviates from the experiment in the upper 
portion of the span and shows good agreement elsewhere. 

Static pressure was measured on the case over the rotor and on the island splitter 
surfaces. A comparison of this experimental data to the numerical solution is shown 
in Fig. 4.15 and Fig. 4.16. For both the island splitter and the case, agreement is 
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Figure 4.12: E cubed bypass stream total pressure ratio 
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Figure 4.15: E cubed circumferential average static pressure: Casing 

seen between the two sets of data. Only on the lower surface of the island splitter is 
a discrepancy shown, and this is probably due again to the mismatched core stream 
performance. 

Static pressure contours for the numerical solution of the fan section are shown in 
Fig. 4.17. Isobar lines are portrayed near the suction and pressure surface with high 
pressure colored red and low pressure colored blue. The shock in the upper portion of 
the rotor is seen in the contour lines transitioning from green to red on each surface. 
The effect the part span shroud has on the flowfield is clearly seen in the figure. 
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Figure 4.16: E cubed circumferential average static pressure: Island splitter 

4.3 Navier Stokes Analysis of Flow Over a Flat Plate 

To verify the boundary conditions and viscous terms modeled in the flow solver 
a laminar boundary layer was simulated. The simulated geometry was a cylindrical 
annulus with a large hub-to-tip radius ratio to minimize curvature effects. A zero 
thickness blade was aligned with the incoming flow, and the blade spacing set so that 
along the local blade surface at midspan, flat plate conditions were well approximated 
(i.e. the blades are far enough apart that there is little effect of the neighboring blade 
on the local boundary layer and midpitch can be thought of as a free stream or infinity 
condition). The grid had 73 axial points and 31 tangential points which were spaced 
so that at midchord of the zero thickness blade, ten points fell in the boundary layer. 
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Figure 4.17: E cubed static pressure contours on solid surfaces 






A laminar boundary layer solution was generated on the grid. The solution was 
compared to a Blasius solution for a boundary layer on a flat plate. A plot of the 
velocity distribution is shown in Fig. 4.18, and the two solutions compare fairly well. 
One region of disagreement is seen about 77 equaling 4, yet as the number of grid points 
in the boundary layer increases (from 9 to 12), the difference between experiment and 
numerical solution decreases. 

4.4 GMA 3007 Fan Section - Navier Stokes Analysis 

While most of the viscous flow solving capability of ADPAC-APES was taken 
from the original code [2] which was validated for turbomachinery flows, a verification 
of the multiple splitter modification was conducted. A Navier-Stokes analysis was 
done on the geometry of the GMA 3007 fan section. The geometry was the same 
as that of the earlier Euler analysis, which was the rotor and the core and bypass 
duct vanes. An axisymmetric plane of the computational grid is seen in Fig. 4.19. 
The full grid for the fan rotor is 141 streamwise points, 59 spanwise points, and 29 
tangential points. The rotor grid models the tip clearance with seven points in the 
gap region. Both the core vane and bypass vane grids are (141X59X23). Table 4.3 
shows pertinent statistics for the grids used. 

The test data for comparison with the numerical solution is the same as that 
shown in the Euler analysis section and the reader is referred to Fig. 4.3 and that 
section for any details. The Navier-Stokes solution obtained was near the design point 
of the fan and was 1.4% higher in flow than the experiment (an improvement from 
3.6% in the Euler solution). 
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Figure 4.19: Axisymmetric plane of GMA 3007 Navier-Stokes grid 




Table 4.3: Navier-Stokes grid parameters for the GMA 3007 geometry 


Blade Row 
Grid 

Streamwise 

Points 

Span wise 
Points 

Tangential 

Points 

Blade Pts 
LE to TE 

Blade Pts 
Hub to Tip 

Fan Rotor 

141 

59 

29 

35 

53 

Core Duct 

141 

59 

23 

27 

19 

Vane 

Bypass Duct 

141 

59 

23 

29 

39 


Vane 


Circumferential- averaged radial profiles of total pressure and temperature ratio 
are shown for the solution and the experimental data in Fig. 4.20. The two sets 
of data show good agreement at this fan condition. The solution compares much 
better to the experimental data than the Euler solution and tnis is probably due to 
several reasons. First, the boundary layer and other viscous phenomena are modeled 
in the Navier-Stokes solution. Also, the resolution of the grid is greater, and the 
tip clearance flow is modeled. Figure 4.21 compares efficiency from the experimental 
data and the numerical solution. Again, the viscous solution compares much better 
with the experimental data than the Euler analysis of Fig. 4.6. 

An effect the tip clearance flow has on the numerical flowfield solution is shown 
in Fig. 4.22. The figure shows the tip region in the fan rotor flowfield. Contours 
of total pressure are shown on an axial plane just downstream of the rotor trailing 
edge. Red contours are the highest total pressure and blue contours are the lowest. 
Streamlines originating from the clearance region just above the blade tip are shown 
as black lines in the figure. Much of the tip clearance flow rolls up into a vortex and 
passes through the contour plane at a low total pressure region. The clearance flow 
clearly effects the local flowfield and its energy makeup. 
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Figure 4.22: GMA 3007 Navier-Stokes solution: Tip clearance flow effects 


5. CONCLUSIONS 


A three-dimensional Euler/Navier-Stokes aerodynamic analysis has been devel- 
oped (through modifying an existing analysis) for fan sections containing multiple ra- 
dial flow splitters. The analysis is capable of calculating internal or external flows and 
can model engine core conditions. The analysis code was verified through comparisons 
with experimental data for two advanced fan section geometries. The numerical solu- 
tions demonstrated reasonable agreement with the experimental data and predicted 
boundary layer characteristics. It was apparent from the comparison to experiment 
that the Navier-Stokes solutions predicted the experimental data better than com- 
parable Euler results. While grid resolution and tip clearance flow are important to 
predictive accuracy, boundary layer modeling has perhaps a stronger effect on predict- 
ing transonic fan performance. The accuracy of the analysis can also be effected by 
additional factors, including errors introduced by turbulence modeling and artificial 
dissipation. 
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